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Abstract 

In  this  research,  an  improved  binary  material-classification  algorithm  is  devel¬ 
oped  to  discriminate  between  metals  and  dielectrics  using  passive  polarimetric  im¬ 
agery  degraded  by  atmospheric  turbulence.  The  technique  implements  a  modified 
version  of  an  existing  polarimetric  blind-deconvolution  algorithm  in  order  to  remove 
atmospheric  distortion  and  correctly  classify  the  unknown  object.  The  classification 
decision  is  based  on  degree  of  linear  polarization  (DoLP)  estimates  provided  by  the 
blind-deconvolution  algorithm  augmented  with  two  DoLP  priors — one  statistically 
modeling  the  polarization  behavior  of  metals  and  the  other  statistically  modeling  the 
polarization  behavior  of  dielectrics.  The  proposed  algorithm  significantly  improves 
upon  a  similar  published  polarimetric  classification  method  by  adaptively  updating 
the  DoLP  priors  as  more  information  becomes  available  about  the  scene. 

Three  approaches  for  implementing  the  adaptive  DoLP  priors  are  presented — 
the  higher-order  super- Gaussian  method,  the  Gaussian  method,  and  the  distribution- 
averaging  method.  The  higher-order  super-Gaussian  method  fits  the  distribution  of 
the  in-progress  DoLP  estimates  from  the  blind-deconvolution  algorithm  with  a  sum 
of  two  super-Gaussian  functions.  The  results  of  the  nonlinear  fit  are  then  used  to 
form  the  DoLP  priors.  The  Gaussian  method  fits  the  distribution  of  DoLP  estimates 
with  the  sum  of  two  Gaussian  functions  to  compute  a  classification  threshold  value. 
The  resulting  threshold  is  then  used  to  update  the  DoLP  priors.  The  distribution- 
averaging  method  approximates  the  threshold  value  by  finding  the  mean  of  the  DoLP 
distribution.  Using  this  threshold  value,  the  DoLP  priors  are  then  formed.  The 
proposed  technique  is  experimentally  validated  by  comparing  classification  results  of 
a  dielectric  and  metallic  sample  obtained  using  the  new  method  to  those  obtained 
using  the  existing  approach.  The  experimental  results  confirm  that  the  new  adap¬ 
tive  method  significantly  extends  the  range  of  validity  of  the  existing  polarimetric 
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classification  technique  to  near-normal  collection  geometries  where  most  polarimetric 
material  classifiers  perform  poorly. 
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Binary  Classification  of  an  Unknown  Object 


through  Atmospheric  Turbulence  Using  a  Polarimetric 
Blind-Deconvolution  Algorithm  Augmented 
with  Adaptive  Degree  of  Linear  Polarization  Priors 

I.  Introduction 

Polarization  is  one  of  the  fundamental  properties  of  light  along  with  irradiance, 
wavelength,  and  coherence  [6].  When  scattered  light  is  reflected  from  an  object, 
its  polarization  signature  provides  information  about  surface  features,  shape,  shading, 
and  roughness  [20,26].  Since  human  activities  often  change  material  surface  features, 
thus  affecting  the  polarization  signature  of  reflected  light,  polarimetric  imaging  has 
emerged  as  a  powerful  tool  to  enhance  understanding  of  an  underlying  scene  of  inter¬ 
est.  For  example,  using  polarimetric  imagery,  Wolff  presented  a  material-classification 
model  to  accurately  discriminate  between  a  dielectric  and  a  metal  based  on  the  po¬ 
larimetric  differences  in  reflecting  light  at  their  surfaces  [27].  This  relatively  new 
polarimetric  imaging  held  has  a  variety  of  applications  ranging  from  remote  sensing 
to  industrial  machine  vision  systems  [26,27]. 

Remote  sensing  refers  to  the  process  of  acquiring  information  about  an  object 
without  physical  contact  with  the  object.  Remote  sensing  can  be  categorized  as  active 
or  passive  depending  on  the  nature  of  the  illumination  source.  In  passive  remote 
sensing,  the  source  of  illumination  is  not  controlled,  e.g.,  the  sun.  An  application 
example  of  passive  remote  sensing  is  the  detection  of  hazardous  industrial  or  warfare 
chemical  agents  [2] .  On  the  other  hand,  active  remote  sensing  uses  a  controlled  source 
to  illuminate  the  object  of  interest.  Ranging  systems  such  as  radio  detection  and 
ranging  (RADAR)  and  laser  detection  and  ranging  (LADAR)  are  common  examples 
of  active  remote  sensing  systems. 
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The  use  of  polarimetric  imagery  has  increased  the  capabilities  of  passive  remote 
sensing  systems,  as  polarization  reveals  information  largely  uncorrelated  with  spectral 
and  irradiance  images  [22,26].  For  example,  Goudail  et  al.  presented  polarimetric 
image  processing  techniques  which  show  that  polarimetric  images  can  improve  target 
detection  and  segmentation  in  remote  sensing  [9].  In  addition  to  Wolff’s  early  work 
in  polarization-based  material  classification,  Thilak  et  al.  introduced  a  polarization- 
based  approach  to  estimate  the  complex  index  of  refraction  for  characterizing  the 
materials  comprising  an  object  [22],  Other  studies  in  material  classification  are  found 
in  [24,28], 

Recently,  Hyde  developed  a  variant  of  the  polarimetric  blind-deconvolution  al¬ 
gorithm  presented  in  [15]  for  classifying  the  material  of  an  unknown  object  using  po¬ 
larimetric  imagery  degraded  by  atmospheric  turbulence  [11].  In  his  algorithm,  prior 
knowledge  about  the  sought  degree  of  linear  polarization  (DoLP),  called  the  DoLP 
prior,  serves  to  classify  an  object  as  a  metal  or  a  dielectric.  Hyde’s  algorithm,  how¬ 
ever,  suffers  when  classifying  an  object  at  near-normal  collection  geometries  because 
of  the  weak  manner  in  which  natural  materials  polarize  light  in  these  conditions.  The 
purpose  of  this  thesis  is  to  present  and  develop  methods  to  alleviate  this  limitation. 

1 . 1  Goals 

The  purpose  of  this  research  is  to  develop  an  algorithm,  somewhat  invariant  with 
respect  to  collection  geometry,  to  discriminate  between  metals  and  dielectrics  using 
polarimetric  imagery  degraded  by  atmospheric  turbulence.  The  quality  of  polarimet¬ 
ric  information  collected  by  a  camera  highly  relics  on  geometry.  The  proposed  method 
enhances  the  poor  performance  of  the  previously  developed  material-classification  al¬ 
gorithm  for  near-normal  collection  geometries  by  adaptively  updating  the  DoLP  priors 
of  a  metal  and  a  dielectric  as  more  information  about  the  target  scene  becomes  avail¬ 
able.  The  method  first  recovers  the  target  scene  from  a  set  of  degraded  polarimetric 
images  using  the  blind-deconvolution  technique.  Applying  the  updated  DoLP  priors, 
the  method  then  discriminates  between  metals  and  dielectrics. 
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1.2  Assumptions 

In  order  to  implement  the  proposed  material-classification  algorithm,  the  fol¬ 
lowing  assumptions  are  made  about  the  passive  polarimetric  remote  sensing  scenarios 
modeled  in  this  research. 

First,  the  unknown  object  in  the  target  scene  is  assumed  to  be  illuminated  with 
a  randomly  polarized  natural  light  source,  e.g.,  the  sun.  Several  studies  have  shown 
that  the  circular  polarization  component  of  light  reflected  from  an  object  illuminated 
with  a  natural  light  source  is  negligible  [4,16,23].  Thus  only  the  linear  components 
of  polarization,  i.e.,  the  Stokes  parameters  So,  Si,  and  S2,  are  considered  throughout 
this  research. 

Secondly,  polarimetric  imagery  is  assumed  to  be  collected  in  the  specular  plane. 
The  specular  plane  is  determined  by  the  directions  of  the  incident  and  reflected  light, 
following  the  simple  geometry  of  reflection  from  a  smooth  (or  mirror-like)  surface  [27]. 
In  passive  polarimetry,  the  angle  of  polarization  (AoP)  provides  information  about 
the  object  geometry,  and  observation  made  in  the  specular  plane  implies  that  the 
AoP  is  zero  [11].  As  a  result,  the  AoP  is  not  considered  in  this  research,  but  should 
be  considered  in  future  research  on  this  topic. 

The  final  assumption  is  that  the  atmospheric  turbulence  is  relatively  weak.  In 
the  observation  plane,  weak  atmospheric  turbulence  randomly  changes  the  light’s 
phase  but  not  its  irradiance  [7].  Strong  atmospheric  turbulence,  on  the  other  hand, 
causes  phase  and  irradiance  fluctuations,  commonly  called  scintillation  [1].  The  algo¬ 
rithm  presented  in  this  research  is  designed  to  correct  phase  aberrations  of  degraded 
images  and  therefore  is  only  applicable  to  weak  atmospheric  turbulence.  This  is  not 
a  very  restrictive  assumption  considering  that  the  light  in  this  analysis  is  assumed  to 
be  spatially  incoherent  (see  assumption  1). 
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1.3  Thesis  Outline 


This  thesis  consists  of  five  chapters  including  the  introduction.  Chapter  II  pro¬ 
vides  the  theoretical  background  necessary  to  understand  the  development  of  the 
material-classification  algorithm.  The  chapter  starts  with  an  explanation  of  polar¬ 
ization  and  imaging  through  atmospheric  turbulence.  A  description  of  the  blind- 
deconvolution  method  is  then  followed  by  a  review  of  the  material-classification  al¬ 
gorithm  in  its  current  form.  Chapter  III  provides  detailed  descriptions  of  the  novel 
methods  for  adaptively  updating  the  DoLP  priors  in  order  to  alleviate  the  afore¬ 
mentioned  limitation  of  collection  geometry.  Chapter  IV  presents  the  results  of  the 
proposed  method  using  polarimetric  imagery  obtained  at  a  near-normal  collection 
geometry.  Finally,  Chapter  V  summarizes  the  findings  of  this  research  and  presents 
possible  topics  for  future  studies. 
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II.  Background 


This  chapter  provides  the  background  necessary  for  understanding  the  polari- 
metric  blind-deconvolution  algorithm  used  for  material  classification.  A  brief 
overview  of  polarization,  incoherent  imaging,  atmospheric  turbulence,  and  blind  de- 
convolution  is  presented. 

2. 1  Polarization 

2.1.1  Polarization  State.  Polarization  simply  describes  the  orientation  of  the 
electric  field.  The  electric  field  E(,s,  t )  at  any  location  z  and  time  t  can  be  represented 
by  two  orthogonal  vector  components,  i.e., 

E  (z,t)  =  E  x(z,t)  +  E  y(z,t)  (2.1) 


where 


E  x(z,  t )  =  E0x  cos  (tat  -  kz) 

E y(z,  t )  =  E0y  cos(c ot  —  kz  —e) 

Here,  E0x  and  E0y  are  the  x —  and  y— component  amplitudes,  u  is  the  angular  fre¬ 
quency,  k  is  the  wavenumber,  and  e  is  the  relative  phase  difference  between  the  vector 
components  [10,20].  Equation  (2.2)  leads  to 
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which  describes  an  ellipse  rotated  through  an  angle  ^  as  shown  in  Fig.  2.1  [10].  The 
angle  ip  is  called  the  AoP  and  is  expressed  as  [20] 


tan 2-0=  cose.  (2.4) 

When  the  relative  phase  difference  e  =  ±7t/2,  ±37t/2,  ±57t/2,  . . . ,  Eq.  (2.3)  reduces  to 
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Figure  2.1:  Illustration  of  the  polarization  ellipse  rotated  through  angle  i/j. 
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which  describes  an  ellipse  in  which  the  major  and  minor  axes  are  aligned  with  the  x 
and  y  axes.  The  resultant  electric  field  vector  E (z,t)  will  appear  to  rotate  and  change 
its  amplitude  as  it  propagates.  Light  in  this  state  is  called  elliptically  polarized  light. 
Elliptical  polarization  is  the  most  general  polarization  state  of  light. 

In  the  case  that  the  x—  and  y— component  amplitudes  are  equal,  i.e.,  E0x  = 
E0y  =  Eq,  Eq.  (2.5)  further  reduces  to 


(E,)2  +  (EJ2  =  El  (2.6) 

yielding  circularly  polarized  light  [10].  If  the  relative  phase  difference  £  —  0,  ±7r,  ±27t,  . . . , 
Eq.  (2.3)  becomes 

E  y  =  ±| °^E0x,  (2.7) 

&0x 
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which  is  the  equation  of  a  line  with  slope  ±E0y/E0x.  Light  in  this  state  is  called 
linearly  polarized  light  [10]. 

2.1.2  Stokes  Vectors  and  Mueller  Matrices.  One  problem  with  modeling 
polarization  using  the  polarization  ellipse  formulation  is  that  it  is  limited  to  scenarios 
involving  fully  polarized  light.  Stokes  vectors  overcome  this  limitation  by  describing 
polarization  in  terms  of  observable  quantities,  namely,  optical  power  or  irradiance  [20]. 
By  taking  time  averages  of  the  electric  held  vector  components,  Eq.  (2.3)  can  be  shown 
to  yield 


where 


C 2  _  c 2  I  q2  i  q2 

°0  —  °1  '  °2  ' 


(2.8) 


So  =  K,  +  el 

c  _  z?2 

-  ^Ox  ~  ^0  y 

52  =  2E0xE0y  COS  £ 

53  —  2E0xE0y  sin  e 


(2.9) 


The  parameter  Sq  is  the  total  irradiance.  The  parameter  S ±  is  the  difference  in  the 
irradiance  of  light  passed  by  a  horizontal  linear  polarizer  and  a  vertical  linear  polarizer. 
If  S'i  =  0,  light  may  be  linearly  or  elliptically  polarized  at  ±45°,  circularly  polarized, 
or  unpolarized  [10].  The  parameter  S2  is  the  difference  in  the  irradiance  of  light  passed 
by  a  linear  polarizer  at  +45°  and  a  linear  polarizer  at  —45°.  Finally,  the  parameter  S3 
is  the  difference  in  the  irradiance  of  light  passed  by  a  right-handed  circular  polarizer 
and  a  left-handed  circular  polarizer.  The  Stokes  parameters  are  often  represented  in 
a  normalized  column  vector  as 


S 


1  S2/S0 


(2.10) 
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The  Mueller  matrix  is  a  4  x  4  matrix  which  models  the  polarization  effects  of 
an  optical  interaction  [20] .  It  relates  the  Stokes  vector  incident  on  an  optical  element 
Sm  to  the  Stokes  vector  leaving  the  element  Sout  as 


Cout 

gout 

s°2ut 

.  S3Ut  . 

Sout  =  MSin 
m00  rn0i  m02  m03 
mw  mu  m12  mi3 
m20  m2i  m22  m23 
m30  m3i  m32  m33 


1 - 

.So 

^0 

st 

.  S a"  . 

(2.11) 


where  M  is  the  Mueller  matrix  of  the  optical  element  [20].  For  example,  the  polar¬ 
ization  effect  on  light  passing  through  an  ideal  linear  polarizer  oriented  at  an  angle  d 
from  the  horizontal  is  expressed  as 


Sout  =  M@Sin 


1 

Co 

oo 

g_ 

1 

cos  29 

sin  2d 

0 

1 - 

So 

Co 

^lout 

l 

cos  2d 

cos22 9 

sin  2d  cos  2d 

0 

^2out 

“  2 

sin  29 

sin  29  cos  29 

sin22d 

0 

st 

1 

O  CO 

CO 

_ i 

0 

0 

0 

0 

\ 

.5  CO 

CO 

_ 1 

S'q11  +  S']11  cos  29  +  St  sin  29 


1  Sq1  cos  29  +  S™cos229  +  S™  sin  29  cos  29 

2  St  sin  29  +  S]11  sin  29  cos  29  +  sin2  20 

0 


(2.12) 


Since  incident  light  is  hltered  through  a  linear  polarizer,  Eq.  (2.12)  shows  that  the  S:i 
component  of  the  output  Stokes  vector  (the  term  describing  the  amount  of  circularly 
polarized  light)  is  zero.  The  total  irradiance  passed  by  the  linear  polarizer  is  contained 
in  the  S0  component  of  the  output  Stokes  vector.  Therefore,  the  irradiance  1(9) 
measured  by  a  detector  with  an  ideal  polarizer  oriented  at  an  angle  9  is  given  by 


1(9)  =  S°ut 


(2.13) 


=  \  (S™  +  Sjn  cos  29  +  S ln  sin  29) . 

Equation  (2.13)  is  an  important  result  and  used  in  the  polarimetric  blind-deconvolution 
algorithm  discussed  in  Section  2.5. 

As  mentioned  in  Chapter  1,  it  is  assumed  that  the  contribution  from  circular 
polarization  is  not  significant  in  this  research.  Therefore,  the  full  Stokes  vector  can  be 
simplified  to  the  linear-polarization-state  Stokes  vector,  eliminating  the  parameter  S% 
whose  value  is  approximately  zero  for  the  cases  of  interest  here.  The  AoP  h  introduced 
in  Eq.  (2.4)  can  be  expressed  in  terms  of  the  Stokes  parameters  as  [20] 

tan2^  =  -2 - -j-  cos  e  —  (2.14) 

A Ox  ~  ^0 y  '-’l 

Like  the  AoP,  the  DoLP  can  be  obtained  from  knowledge  of  the  Stokes  parameters 
as  well: 


DoLP  Jj*  =  dgpl  (215) 

hot  *So 

where  Ipoi  is  the  linearly  polarized  light  irradiance  and  hot  is  the  total  irradiance.  Its 
range  is  between  0,  i.e.,  randomly  or  unpolarized  state,  and  1,  fully  polarized. 


2. 2  Imaging 

2.2.1  Incoherent  Imaging.  Incoherent  illumination  implies  that  each  spatial 
point  on  the  object  field  is  statistically  independent  of  all  other  points  [18].  In  the 
case  of  incoherent  imaging,  the  image  irradiance  i(x,y)  can  be  found  by  convolving 
the  object  irradiance  o(x,  y)  and  the  irradiance  impulse  response  of  the  optical  path 
of  interest,  also  known  as  the  point  spread  function  (PSF)  \h(x,  y)\2: 


i(x,y)  = 


|  h(x  —  u,y  —  v)\2o(u,v)dudv 


(2.16) 
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where  (x,y)  defines  an  image  plane  coordinate  and  (u,v)  defines  an  object/target 
plane  coordinate.  Note  that  because  of  the  nature  of  imaging  (u,  v )  and  ( x ,  y)  are 
typically  related  by  the  magnification  M.  Without  loss  of  generality,  unit  magnifica¬ 
tion  is  assumed  in  all  analysis  to  follow.  Taking  the  Fourier  transform  of  Eq.  (2.16) 
gives  the  equivalent  relationship  of  incoherent  imaging  in  the  spatial-frequency  domain 
as 


I(/»,/y)  =  H(/x,/w)0(/„/y)  (2.17) 

where  I(fx,  fy),  0(/x,  fy),  and  H (fx,  fy)  are  the  Fourier  transforms  of  i(x,y),o(x,y), 
and  | h(x,  y)  |2,  respectively,  and  (fx,  fy)  is  the  spatial-frequency  coordinate  correspond¬ 
ing  to  the  image-plane  coordinate  (x,  y )  [8].  Normalizing  H (fx,  fy)  by  its  value  at  zero 
frequency  yields  the  function  H  [8]: 


n(fx,fy) 


I IZ  \h  Of  y)  I2  exP [-jZnjfxX  +  fyy)]dxdy 
I  IZo  \h(x,y)\2dxdy 


(2.18) 


The  function  'H  is  called  the  optical  transfer  function  (OTF),  and  its  modulus  \H\  is 
referred  to  as  the  modulation  transfer  function  (MTF),  commonly  used  to  characterize 
imaging  systems. 


2.2.2  The  OTF  of  a  Diffraction- Limited  Incoherent  Imaging  System.  An 
optical  system  is  said  to  be  diffraction  limited  if  an  image  is  limited  in  quality  only  by 
diffraction,  not  by  imperfections  in  the  imaging  system  [8].  The  OTF  of  a  diffraction- 
limited  incoherent  imaging  system  with  a  circular  aperture  of  diameter  D  is  given 
by  [8] 


2 


n{P)  =  { 71 


o 


if  P<  2p0 

1 


otherwise 


(2.19) 
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where  p0  is  the  cutoff  frequency  of  a  coherent  imaging  system  of  focal  length  /  and 
wavelength  A,  i.e., 


Po  =  (2.20) 

The  range  of  the  radial  distance  p  =  \J +  fy  in  Eq.  (2.19)  indicates  that  an  incoher¬ 
ent  imaging  system  extends  the  cutoff  frequency  to  twice  that  of  a  coherent  imaging 
system  [8]. 

2.2.3  Atmospheric  Turbulence.  Atmospheric  turbulence  is  predominantly 
caused  by  temperature  fluctuations  in  the  atmosphere  resulting  from  differential  heat¬ 
ing  and  cooling  of  Earth’s  surface  by  the  sun  [18].  These  temperature  fluctuations 
cause  spatial  and  temporal  fluctuations  in  the  index  of  refraction  of  the  atmosphere 
resulting  in  random  spatial  and  temporal  variations  in  optical  path  length  (OPL)  ex¬ 
perienced  by  transiting  light  [18].  These  OPL  variations  manifest  phase  aberrations 
in  the  optical  wavefront  as  the  light  propagates  through  the  atmosphere.  These  phase 
aberrations  add  another  layer  of  distortion  in  addition  to  the  imaging  system.  As  a 
result,  the  quality  of  received  images  is  degraded  and  the  resolution  of  the  imaging 
system  is  reduced. 

2.2.4  Imaging  through  Turbulence.  The  instantaneous  PSF  and  its  OTF 
through  atmospheric  turbulence  are  random  due  to  the  stochastic  nature  of  the  atmo¬ 
sphere,  as  described  above  [18].  Physical  meaning  can  only  be  gleaned  by  inspection 
of  the  average  PSF  and  its  Fourier  transform,  the  average  OTF.  To  see  how  atmo¬ 
spheric  phase  aberrations  are  modeled  using  the  OTF  TL,  Eq.  (2.18)  is  first  expressed 
in  a  different  form  [7,8]: 


n(fx,fy) 


f  IZ  H(x>  y)H*(x  -  v  -  ^ffy)dxdv 

f  fZo  \H(x,y)\2dxdy 


Here,  the  function  H  represents  a  modified  aperture  function,  i.e., 


(2.21) 
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H{x,y)  =  P(x,y)t(x,y) 


(2.22) 


where  P(x,y )  is  the  classic  aperture  apodization  function  and  t(x,y )  is  the  complex 
transmittance  given  by 


t(x,y)  =  exp\jd(x,y)\.  (2.23) 

In  Eq.  (2.23),  9(x,y )  represents  the  random  phase  induced  by  the  atmosphere  [7]. 
Substituting  Eq.  (2.22)  into  Eq.  (2.21)  yields  an  expression  for  the  instantaneous 
OTF  of  an  imaging  system  including  atmospheric  effects  7-Lsys  [7]: 


'Hsys(fx,fy 


f  JZcP(X’y)P*(X  -  Xffx,y~  Xffy)t(x,y)t*(x 

J  JZJp(x^y)\2dxdy 


A ffx,y~  A ffy)dxdy 


(2.24) 


The  average  OTF  of  the  imaging  system  ' Hsys{fx ,  fy)  is  computed  by  taking  the 
expectation  of  Eq.  (2.24)  yielding 


' HsysifxJy )  = 


f  fZoP(x’y)P*(x  -  Xffx,y-  A ffy)E  [t(x,y)t*(x  -  Xf  fx,y  ~  Xf  fy)]  dxdy 


HZ  I p(x,y)\2dxdy 


(2.25) 


Equation  (2.25)  can  be  written  as 


nsys(fx,  fy)  =  Ho  Ur.  fy)  Us  (/,•  fy) ,  (2.26) 

where  ' Ho(fx ,  fy)  is  the  OTF  of  a  diffraction-limited  system  and  Us(fx,  fy)  represents 
the  average  OTF  of  the  induced  random  phase  [7].  Two  types  of  the  average  OTF 
are  discussed  in  this  section — long-exposure  and  short-exposure  OTFs. 
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Long-exposure  imaging  is  imaging  performed  over  a  long  integration  time.  Be¬ 
cause  of  the  long  integration  time,  many  instances/realizations  of  atmospheric  turbu¬ 
lence  are  captured  in  the  resulting  image.  The  long-exposure  OTF  of  the  atmosphere 
was  derived  by  Fried  as 


1-Lle  (p)  =  exp 


(2.27) 


where  r0  is  the  atmospheric  coherence  diameter  [19].  The  atmospheric  coherence 
diameter  describes  the  largest  effective  aperture  diameter  for  image  resolution.  For  a 
plane  wave  source,  r0  can  be  computed  as 


0.423 k2 


-3/5 


(2.28) 


where  C\  is  the  refractive-index  structure  parameter,  which  is  a  measure  of  the  tur¬ 
bulence  strength,  and  L  is  the  propagation  distance  [19]. 


Short-exposure  imaging  is  imaging  performed  over  a  short  integration  time.  In 
short-exposure  imaging,  only  a  single  instance  of  atmospheric  turbulence  is  captured 
in  the  image  thus  effectively  freezing  the  effects  of  the  atmosphere  [18].  The  key  differ¬ 
ence  between  short-exposure  and  long-exposure  imaging  is  that  images  taken  using  a 
short  exposure  are  not  distorted  by  wavefront  tilt  which  has  the  effect  of  broadening 
the  PSF  (or  equivalently  narrowing  the  OTF)  in  the  long-exposure  case  [18].  The 
short-exposure  OTF  of  the  atmosphere  is  given  by 


T-isE  (p)  =  exp 


a 


(2.29) 


where  a  is  1/2  or  1  when  scintillation  is  or  is  not  present,  respectively  [19].  Note 
that  Use  reduces  to  Hle  by  setting  a  =  0.  Figure  2.2  illustrates  how  atmospheric 
turbulence  corrupts  an  image  and  reduces  the  resolution  of  an  imaging  system  as 
compared  to  the  diffraction-limited  case. 
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(a)  True  object 


(b)  Instantaneous  short-exposure  PSF 


(c)  Long-exposure  PSF 


# 


(d)  Diffraction-limited  image 


(e)  Instantaneous  short-exposure  image 


(f)  Long-exposure  image 


Figure  2.2:  Effect  of  atmospheric  turbulence  on  imaging. 

2.3  Blind  Deconvolution 

Blind  deconvolution  refers  to  the  process  of  estimating  the  true  object  and 
the  PSF  from  degraded  images  with  no  or  partial  information  about  the  imaging 
system  or  object  [14].  Blind  deconvolution  has  applications  in  various  fields  including 
optics,  astronomy,  signal  processing,  and  seismology.  There  exist  many  methods  to 
solve  the  blind-deconvolution  problem,  such  as  zeros  sheet  separation,  a  priori  blur 
identification,  parametric  estimation,  etc.  [14].  Among  these  approaches,  the  most 
relevant  to  this  thesis  is  the  maximum-likelihood  (ML)  estimation  method,  which  is 
a  parametric  estimation  technique.  ML  blind-deconvolution  algorithms  developed  by 
Schultz  and  LeMaster  and  Cain  will  be  summarized  here  after  a  brief  review  of  ML 
estimation. 
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2.3.1  Maximum- Likelihood  Estimation.  The  ML  method  solves  a  blind- 
deconvolution  problem  by  estimating  the  parameter  set  0  =  [di,d2, . . .  ,6n]t  of  a 
distribution  X  based  on  the  observed  data  x  =  [aq,  x2, . . . ,  xn]T  drawn  according 
to  the  distribution  X  [17].  An  estimate  of  the  parameters  is  made  to  maximize  the 
probability  or  likelihood  of  observing  data  x  given  the  parameter  set  6  [14].  Let 
/(x|0)  be  the  conditional  probability  density  function  (PDF)  of  the  distribution  X 
given  the  parameter  set  6  [17].  The  ML  estimate  of  0  is  then 

0rni  =  arg  |  max  /x  (»)}  =  arg  |  max  /(x  l«)},  (2.30) 

where  lx(0)  denotes  the  likelihood  function  [14].  Since  most  distributions  of  engineer¬ 
ing  interest  belong  to  the  exponential  family  such  as  Gaussian,  Rayleigh,  Poisson,  etc., 
the  logarithm  of  the  likelihood  function  is  often  utilized  in  Eq.  (2.30)  [17].  Note  that 
the  logarithm  is  a  monotone  function,  thus  maximizing  the  log-likelihood  function  is 
equivalent  to  maximizing  the  likelihood  [17]: 

L(0)  =  \oglx(0)  =  log/(x|0).  (2.31) 

2.3.2  Expectation- Maximization  Algorithm.  The  expectation-maximization 
(EM)  algorithm  is  an  iterative  algorithm  to  produce  ML  estimates  of  the  parameter 
set  6.  In  the  EM  algorithm,  the  observed/measured  data  y  =  [yi,  y2, . . . ,  |/at]T,  termed 
the  incomplete  data,  is  mapped  to  a  new  set  of  data  x  =  [aq,  ;r2, . . . ,  xn}t ,  termed  the 
complete  data,  where  the  ML  estimation  problem  is  more  easily  solved  [21].  Note  that 
the  complete  data  x  are  not  observed  directly  and  thus  do  not  need  to  have  a  physical 
meaning.  They  only  need  to  be  statistically  consistent  with  the  incomplete  data  y  [15]. 
The  idea  behind  the  EM  algorithm  is  that  even  though  the  complete  data  specification 
is  not  known,  making  use  of  the  underlying  distribution  log[/(x|0)],  an  estimate  of 
6  can  be  determined  [17].  Obtaining  this  estimate  Qmi  is  a  nonlinear  optimization 
problem  of  several  variables  [14].  Among  various  implementation  methods,  the  EM 
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algorithm  is  the  most  popular  for  its  simplicity  and  computational  efficiency  [14]. 
The  EM  algorithm  converts  the  nonlinear  optimization  problem  into  a  simpler  linear 
iterative  procedure,  making  implementation  straightforward. 

The  goal  of  the  EM  algorithm  is  to  find  the  parameter  set  0  which  maximizes 
the  log-likelihood  log[/(x|0)]  [17].  Recall  that  the  complete  data  x  are  not  directly 
observed  and  thus  it  is  not  possible  to  compute  the  log-likelihood  of  x  directly  [17]. 
Instead,  the  EM  algorithm  maximizes  the  expectation  of  log[/(x|0)]  given  the  incom¬ 
plete  data  y  and  the  current  estimate  of  the  parameter  set  0  [17].  The  expectation 
step  (the  first  step  of  the  EM  algorithm)  computes  this  conditional  expectation,  i.e., 


Q(0|0'i)=£[log/(x|0)|y,e'i]  (2.32) 

where  Q(6\6k )  is  called  the  objective  function  and  6k  represents  the  estimate  of  0 
after  k  iterations  [17].  Once  the  objective  function  Q  is  defined,  the  maximization 
step  (the  second  step  of  the  EM  algorithm)  provides  a  new  estimate  of  the  parameter 
set  0k+1  by  maximizing  Q  [17], 


6k+1  =  arg  {maxQ(0|0fc)}  .  (2. 

The  EM  algorithm  starts  with  choosing  an  initial  0  estimate  and  then  is  iterated  until 
convergence  or  a  stopping  criterion  is  met. 

2.f  Multiframe  Blind-Deconvolution 

Unconstrained  ML  estimation,  discussed  in  Section  2.3,  can  converge  to  a  trivial 
solution,  namely,  the  object  is  estimated  as  a  point  source  and  the  PSF  is  estimated 
as  the  measured  data.  This  estimate  needs  to  be  avoided  because  it  provides  no  infor¬ 
mation  about  the  true  object  [21].  Schulz  introduces  two  techniques  for  performing 
ML  blind  deconvolution  of  a  sequence  of  short-exposure  images  while  avoiding  the 
trivial  solution  [21].  The  first  method,  penalized  ML  estimation,  adds  a  penalty  term 
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to  the  original  log-likelihood  function  so  that  point-source-like  object  estimates  are 
penalized  [21].  The  second  proposed  method,  which  is  applicable  to  this  research, 
performs  ML  blind  deconvolution  using  a  parameterized  form  of  the  PSF  [21].  This 
approach  is  described  in  more  detail  below. 

Let  o(x)  and  hk( y|x)  denote  the  object’s  irradiance  and  the  PSF  for  the  kth 
short-exposure  image,  respectively.  The  coordinate  x  =  {x\,xi)  defines  a  point  in  the 
object/target  plane;  the  coordinate  y  =  (y  1,7/2)  defines  a  point  in  the  image  plane 
[11,21],  Assuming  discretized  object  and  image  plane  regions,  the  ideal  irradiance  in 
the  kth  frame  is 


4(y;  o,  hk)  =  22  /ifc(y|x)o(x).  (2.34) 

X 

This  ideal  irradiance  is  never  actually  measured.  What  is  actually  detected 
is  a  noise-degraded  version  of  ik.  The  assumption  of  this  algorithm  and  the  other 
blind-deconvolution  algorithms  to  follow  is  that  shot  noise  is  the  dominant  noise 
source  degrading  ik.  Let  dk( y)  be  defined  as  the  data  detected  at  location  y  in  the 
kth  frame  where  ik(y,o,  hk)  =  E[dk(y)].  Since  shot  noise  is  Poisson  distributed,  the 
log-likelihood  L(o,  h )  becomes 


L(o,  h)  =  Y Z 

k 


y^4(y;  o,  hk)  +  22  4(y)  ln*fc(y;  o,  hk) 

y  y 


(2.35) 


after  eliminating  terms  that  do  not  affect  the  maximization  [21],  The  PSF  is  param¬ 
eterized  by  phase  errors  over  the  aperture  as 


hk{ y|x;  ak,  9k)  =  akg( y  -  x;  9k), 


(2.36) 


where 
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2 


(2.37) 


<?(x;  4) 


7L(u)  exp  [i0fc(u)]  exp(— *27tku  •  x)  , 

U 


7l(u)  is  the  aperture  function,  u  denotes  a  spatial  coordinate  pair  in  the  pupil  plane, 
ak  is  the  gain  of  the  kth  PSF,  9k  models  the  phase  aberrations  caused  by  atmospheric 
turbulence  or  the  image  system  in  the  kth  frame,  and  k  is  a  scaling  constant  between 
pupil-plane  and  image-plane  coordinates  [21]. 

The  EM  algorithm  is  then  used  to  optimize  estimates  of  the  object  irradiance 
and  the  parameters  ak,9k-  Substitution  of  Eqs.  (2.34)  and  (2.36)  into  Eq.  (2.35),  the 
log  likelihood  becomes 


L(o,  h )  =  —  EEE  akg(y  -  x;  9k)o(x)  +  EE  4(y)  hr  oik 

k  y  X  k  y 

+EE  4(y)  In  P  g(y  -  x;  9k)o(x) 

k  y  Lx 


(2.38) 


The  ML  estimate  of  ak  can  be  obtained  by  differentiating  Eq.  (2.38)  with  respect  to 
ak  and  setting  the  derivative  to  zero  [21].  The  resulting  estimate  ak  is 


ak  =  Dk/G ,  (2.39) 

where  Dk  =  4( y)  and  G  =  Xlxfi,(x>  4)  [21].  The  update  equations  for  the  object 

irradiance  o(x)  and  the  phase  parameter  9k  are  derived  by  applying  the  method  of 
Lagrange  multipliers  and  the  Gerchberg-Saxton  algorithm,  respectively: 


onew( x)  =  zrVw(x)  Y 

k  y 


akg(y  -  x;  9°ld) 
4(y;  oold ,  hlld) 


dk  (y) 


(2.40) 


and 
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[  9{u)  if  Of1)  In  g(x\  9)  >  J^£(x;  Of1)  ln#(x;  Of1) 

er(u)=l  x  x  ,  (2.41) 

[9°kld(u)  else 

where 

6(u)  =  phase  ^F”1  |^(x;  ff)  exp  [iphase  {<?(x;  9°fd)}]  }) 

fyOld  f  -it-  _  -%r  \ 

«x;  Of)  =  9(x;  Of)  £  -  <™0 

<Kx; o°kd )  =  F*  {^(«) exP  [^(u)] } 

Here,  Fz  denotes  a  discrete  Fourier  transform  [21],  Figure  2.3  depicts  the  ML  estima¬ 
tion  results  using  Schulz’s  algorithm  just  described.  The  first  row  of  Fig.  2.3  is  the 
true  object.  The  second  row  shows  the  two  short-exposure  PSFs  used  in  this  simula¬ 
tion  ( D/r0  =  20).  The  third  row  contains  the  detected  images,  which  are  corrupted 
by  atmospheric  turbulence  and  simulated  shot  noise  (signal-to-noise  ratio  (SNR)  = 
\/l000).  The  fourth  row  is  the  estimated  object  after  200  iterations.  The  signifi¬ 
cance  of  the  result  is  that  the  object  estimate  clearly  indicates  two  distinct  objects  of 
different  sizes.  On  the  contrary,  this  information  is  partly  lost  in  the  detected  images. 

2.5  Multichannel  Blind  Deconvolution  of  Polarimetric  Imagery 

Recently,  LeMaster  and  Cain  developed  a  polarimetric  version  of  Schulz’s  polar¬ 
ization  insensitive  algorithm.  In  their  study,  the  linear  polarization  state  is  considered, 
and  the  images,  measured  by  detectors  at  different  polarization  orientations  (chan¬ 
nels),  are  fused  to  provide  polarimetric  information  about  the  scene,  namely  estimates 
of  the  AoP  as  well  as  estimates  of  the  object  and  the  channel  PSFs  [15]. 


19 


Figure  2.3:  ML  estimation  results  using  Schulz’s  parameterized-PSF  blind- 

deconvolution  algorithm.  The  first  row  is  the  true  object.  The  second  row  contains 
the  two  short-exposure  PSFs  used  in  this  simulation  (D/ro  =  20).  The  third  row 
shows  the  detected  images  (SNR  =  \/l000).  The  fourth  row  is  the  estimated  object 
after  200  iterations. 

As  reviewed  in  Section  2.1.2,  the  irradiance  measured  at  a  detector  with  a  linear 
polarizer  orientated  at  angle  9  is  given  by 


m 


1 

2 


(Sq  +  Si  cos  29  +  S2  sin  29) 


(2.43) 
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where 


51  =  PSo  cos  2 a 

52  =  PS(i  sin  2a 


(2.44) 


Here,  a  is  the  AoP  and  P  is  the  DoLP.  Substitution  of  Eq.  (2.44)  into  Eq.  (2.43)  and 
subsequent  simplification  yields 


oc(x)  =  ^A„(x)  +  Ap(x)cos2  (a(x)  -  9C) .  (2.45) 

where  \u  —  (1  —  P)Sq  and  \p  =  PSo  represent  the  unpolarized  and  polarized  compo¬ 
nents  of  the  scene,  respectively  [15].  Given  the  PSF  hc  and  the  object  irradiance  oc, 
the  ideal  irradiance  ic  is 


*c(y)  =  °c(x)/rc(y  -  x).  (2.46) 

X 

As  in  Schulz’s  algorithm,  the  incomplete  data  dc( y)  denote  the  number  of  pho¬ 
tons  at  location  y  in  polarization  channel  c.  Each  dc  is  subsequently  split  into  polar¬ 
ized  and  unpolarized  components  dpc{ y,x)  and  duc( y,  x),  which  in  accordance  with 
the  dictates  of  the  EM  algorithm  become  the  complete  data: 

dc{  y)  =  x)  +  ^(y, x)  (2-47) 

X  X 

where 


1 

E[dpc(  y,x)]  =  -Au(x) 

E[duc{ y,x)]  =  Ap(x)cos2  (a(x)  -  9C) 


(2.48) 


Substituting  Eqs.  (2.45)-(2.47)  into  Eq.  (2.35),  the  complete  data  log-likelihood  func¬ 
tion  Lcd  becomes  [15] 
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c  y  x 


LCD(XU ,  A p,a,/i)  =  XiEE  dwc  (y, x)  In  [|Aw(x)/ic  (y  -  x)]  -  §A tt(x)/ic  (y  -  x) 


EEEi^c  (y,  x)  In  [Ap(x)cos2(a(x)  -  6c)hc  (y  -  x)] 

c  y  x  ^ 

-Ap(x)cos2(a(x)  -  6c)hc  (y  -  x)} 


(2.49) 


Update  equations  for  the  sought  parameters  are  obtained  by  calculating  the 
conditional  expectation  of  LCD  (expectation  step)  and  maximizing  the  resulting  ex¬ 
pression  for  Au(x),  Ap(x),  a(x),  and  /ic(y|x).  The  resulting  update  equations  for  the 
polarization  components  Xu  and  Xp  at  pixel  xq  are  [15] 


<2-50) 

c  y 

where  k  refers  to  either  u  or  p  for  the  unpolarized  or  polarized  component,  respectively, 
C  is  the  number  of  channels,  and 

V£+1(y,x)  =  4rk‘(x)“s2(“"(x)  -  »c)K  (y  -  *) 

\cv)  s  .  (2.51) 

Cc+1(y>x)  =  (y  -  x) 

2  w) 

The  AoP  estimate  is  [15] 


(2'52) 

Finally,  the  PSF  is  estimated  according  to  the  update  equations  introduced  by  Schulz 
[Eqs.  (2.36),  (2.37),  (2.41),  and  (2.42)],  where  the  gain  ak  is  set  to  unity  such  that 

hk{ y|x;  ak,  9k )  =  g( y  -  x;  9k).  (2.53) 
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(a)  Object  (b)  Unpolarized  (c)  Polarized 


Figure  2.4:  ML  estimation  results  using  three-channel  (0°,  60°,  120°)  polarimetric 
imagery.  The  first  row  contains  the  true  object  and  its  unpolarized  and  polarized 
components.  The  second  row  depicts  the  simulated  short-exposure  PSF  (D/tq  =  20). 
Shown  in  the  third  row  are  the  “measured”  S0  (SNR  =  \/l000),  and  the  estimated 
unpolarized  and  polarized  components  of  the  scene  after  200  iterations  of  LeMaster 
and  Cain’s  algorithm 

Figure  2.4  depicts  estimation  results  using  simulated  polarimetric  imagery  col¬ 
lected  using  three  polarimetric  channels.  The  analyzers  in  the  three  channels  are 
oriented  at  0°,60°,  and  120°.  The  target  scene  has  two  objects — one  fully  polarized 
and  the  other  unpolarized,  as  shown  in  the  first  row  of  Fig.  2.4.  These  objects  are 
imaged  through  atmospheric  turbulence  whose  short-exposure  PSF  (D/r o  =  20)  is 
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shown  in  the  second  row.  The  third  row  contains  the  “measured”  S0  with  SNR  = 
\/l000.  Also  shown  in  the  third  row  are  the  unpolarized  and  polarized  components 
of  the  scene  estimated  after  200  iterations.  Note  that  the  object  and  its  polarimetric 
content  are  fully  recovered. 

2. 6  Classifying  the  Material  of  an  Unknown  Object  with  DoLP  Priors 

This  section  reviews  a  material  classification  method  using  degraded  polari¬ 
metric  imagery  introduced  by  Hyde.  Using  a  variant  of  the  polarimetric  ML  blind- 
deconvolution  algorithm  developed  by  LeMaster  and  Cain,  Hyde  estimates  the  true 
object’s  irradiance  So,  the  DoLP,  the  AoP,  and  the  PSF  [11].  Unlike  LeMaster  and 
Cain,  Hyde  adds  prior  knowledge  about  the  sought  DoLP  estimate,  called  the  DoLP 
prior,  into  his  algorithm  to  serve  as  a  means  to  classify  the  unknown  object.  This 
DoLP  prior  distribution  is  simply  referred  to  as  the  prior  hereafter.  Hyde  introduces 
two  priors — one  representing  dielectric  materials  and  the  other  representing  metallic 
materials.  These  priors  are  discussed  in  more  detail  in  the  next  section,  ft  should 
be  noted  that  by  introducing  prior  information  about  the  sought  parameters  into  the 
log-likelihood  function,  the  problem  changes  from  being  an  ML  estimation  problem 
to  a  maximum  a  posteriori  (MAP)  problem.  Since  MAP  estimation  is  very  similar 
in  nature  to  ML  estimation  (previously  reviewed),  MAP  estimation  is  not  reviewed 
here.  The  interested  reader  is  referred  to  Refs.  [25]  and  [3]  for  more  information. 

2.6.1  DOLP  Priors.  Hyde  adds  two  DoLP  priors  into  his  algorithm  to 
classify  dielectric  and  metallic  materials  [11,12],  To  formulate  these  priors,  he  predicts 
the  DoLPs  of  dielectrics  and  metals  using  a  polarimetric  bidirectional  reflectance 
distribution  function  (pBRDF)  [11],  These  pBRDF  results  show  that  dielectric  and 
metal  DoLPs  are  uniformly  distributed  within  certain  nearly  exclusive  ranges  [11].  He 
approximates  the  non-differentiablc  uniform  distribution  using  a  continuous  super- 
Gaussian  distribution  as 
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Figure  2.5:  DoLP  priors  modeled  with  a  super-Gaussian  distribution.  The  con¬ 

stants  in  Eq.  (2.54)  are  c  —  1,  a  —  7.2,0  =  0.15,  m  =  10  for  the  metal  DoLP  prior 
and  c  =  1,  a  =  3.1, 0  =  0.65,  m  —  10  for  the  dielectric  DoLP  prior  [11], 


n(-P)  =  cexp  {  — [a(P  —  0)}m}  ,  (2.54) 

where  P  denotes  the  DoLP,  c  is  a  normalizing  constant  to  ensure  that  IT(P)  integrates 
to  one,  a  is  a  constant  equivalent  to  the  variance  for  the  super-Gaussian  distribution,  0 
is  a  constant  which  controls  the  center  of  the  distribution,  and  m  is  an  even  integer  [11]. 
Figure  2.5  depicts  the  DoLP  priors  used  by  Hyde  in  [11]. 

The  new  complete-data  log-likelihood  function  incorporating  the  DoLP  prior  is 


\(l~P)S0hc  (y-x) 

-\(l~  p)soK  (y  -  x)|  +  {<Wy> x) ln  [^v°s2(«(x)  -  °c)  (2.55) 

'  c  y  x 

xhc  (y  -  x)]  — PS'0cos2(a(x)  -  0c)hc  (y  -  x)}  +  EEEMn(p>i 


L(P,S0,a,h)  = 
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where  duc( y.x)  and  dpc( y.x)  are  the  incomplete  data  as  defined  in  Eq.  (2.47)  [11].  In 
Hyde’s  algorithm,  the  AoP  a(x)  is  assumed  to  be  zero,  which  implies  that  observation 
is  made  in  the  specular  plane  [11]. 

As  in  Schulz’s  and  LeMaster  and  Cain’s  algorithms,  Hyde  employs  the  EM 
algorithm  to  find  estimates  of  the  desired  parameters  P,So,a,  and  hc.  The  outcome 
of  the  expectation  step  is  the  objective  function  Q  which  is  [11]: 


Qn+\p , s0, a, h)  =  Y, E E \E  [<Uy> x 


c  y  X 


In 


-{l  -  Pn)S%h™  (y-x) 


-3(1-  pn)SoK  (y -x)  +  JEE{E  [4(y,  x)J  In  [Pn  SqCos2  (a11  (x.)  -  0e) 

'  c  y  x 

xhnc  (y  -  x)]  —PhSq cos2 (an (x)  -  9c)hnc  (y  -  x)}  +  £££  In  [n(Pn)] 


c  y  X 


(2.56) 


Differentiating  Q  with  respect  to  P,  So,  a,  and  hc  and  setting  the  derivatives  equal  to 
zero  yield  the  update  equations 


0  =  P(1-P)££ 


1  dn(p) 
n(P)  dp 


-p  (Cc+1( y, x)  +  V£c+1(  y, x))  +  XI ^c+1(y, x) 

c  y  c  y 

*»  =  ^££  (Cc+1(y> x)  +  ^pC+1(y> x)) 


c  y 


0  =  ^+1  tan(a(x)  -  6C 


c  y 


Mz)  =  -jj £  [V>yl(y.y-z)  +  V>pC+1(y,y-z 


(2.57) 


where 
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1  dr 


Cc  (y,x)  =  E  duc(y,*)\dc,Pn,S%,an,hnc  =  -^(1  -P»)^ 

-I  z 


VC+1(y,x)  =  E  4,(y,x)K,P",sy,a”,ft; 


=  -P"S?cos2(a" 
u  v 


(2.58) 


z  =  y  —  x,  and  _DC  =  ^2  dc( y)  [11].  The  channel  PSF  hc( z)  is  estimated  according  to 
Schulz’s  parameterized  PSF  technique  introduced  in  Section  2.4.  Applying  the  DoLP 
prior  model  II(P)  in  Eq.  (2.54)  to  the  DoLP  update  equation  in  the  first  row  of  Eq. 
(2.57)  yields  a  DoLP  update  equation  of  the  form  [11] 


0  =  mamP2(P  -  ft)™-1  -  mamP(P  -  ft)™-1 

[«o+1(y.x)+^c+1(y.x)l +E5>?‘(y-x) '  <2,59) 

C  y  C  y 

To  classify  each  pixel,  two  sets  of  estimates  for  the  parameters  P,  So,  and  h  are 
computed  from  Eq.  (2.57),  given  the  metal  and  dielectric  priors  [11].  These  two  sets 
of  estimates  are  then  substituted  into  the  objective  function  Q  which  is  simplified 
in  Eq.  (3.1)  using  Eqs.  (2.56)  and  (2.58).  Whichever  set  of  estimates  maximizes  Q 
classifies  that  pixel  either  as  a  metal  or  as  a  dielectric. 

Figure  2.6  depicts  the  recovery  of  the  target  scene  from  turbulence-degraded 
(. D/r0  =  20)  images  using  four  polarimetric  channels,  i.e.,  9  =  0°,45°,90°,  135°.  The 
first  row  of  the  figure  is  the  true  target  scene,  which  consists  of  an  object  composed  of 
metallic  and  dielectric  components.  The  portions  of  the  object  with  the  lower  DoLP 
P  =  0.2  represent  metallic  components;  the  portions  with  the  higher  DoLP  P  =  0.7 
represents  dielectric  components.  The  “measured”  So  (SNR  =  v/ 1000)  and  DoLP  are 
shown  in  the  second  row.  The  third  row  contains  the  estimates  of  the  scene,  applying 
the  update  equations  in  Eq.  (2.57),  after  200  iterations.  Figure  2.7  shows  the  material 
classification  results,  incorporating  Hyde’s  metal  and  dielectric  DoLP  priors.  As  Figs. 
2.6  and  2.7  indicate,  the  algorithm  successfully  recovers  and  classifies  the  target  scene. 
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(b)  True  DoLP 


(a)  True  SQ 


Figure  2.6:  Estimation  results  using  four-channel  (0°,  45°,  90°,  135°)  polarimetric 

imagery.  The  top  row  is  the  true  target  scene.  The  second  row  shows  measurements 
(SNR  =  VlOOO)  degraded  by  atmospheric  turbulence  {D/tq  =  20).  The  third  row 
contains  the  final  estimation  results  after  200  iterations. 

2.6.2  Limitations.  In  Hyde’s  algorithm,  fixed  DoLP  priors  are  used  to  clas¬ 
sify  the  unknown  object.  A  key  limitation  of  this  approach  lies  in  the  fact  that  natural 
materials  weakly  polarize  scattered  light  at  near-normal  incident  angles  [11].  As  a  re¬ 
sult,  the  images  collected  at  near-normal  geometries  (0j,  6r  =  24°  in  Hyde’s  paper)  are 
poorly  classified  by  his  material-classification  algorithm  [11],  Hyde  suggests  that  this 


Figure  2.7:  Binary  material  classification  of  an  unknown  object  (red  =  metal,  blue 
=  dielectric). 

limitation  can  be  alleviated  by  adaptively  updating  the  DoLP  priors  as  more  infor¬ 
mation  becomes  available  about  the  scene  [11].  Chapter  III  presents  a  methodology 
to  implement  the  adaptive  DoLP  priors  suggested  by  Hyde  and  represents  the  main 
theoretical  thrust  of  this  research. 


29 


III.  Research  Methodology 


This  chapter  describes  a  methodology  to  implement  a  polarization-based  material- 
classification  algorithm  incorporating  adaptive  DoLP  priors.  The  process  in¬ 
volves  updating  the  DoLP  priors  by  inspecting  the  distribution  of  DoLP  estimates  to 
enhance  classification  results  at  near-normal  collection  geometries. 

3. 1  Material  Classification 

The  material-classification  algorithm  introduced  by  Hyde  classifies  the  materi¬ 
als  of  an  unknown  object  using  polarimetric  imagery  degraded  by  the  atmosphere. 
Incorporated  in  his  algorithm  are  DoLP  priors  of  a  metal  and  a  dielectric.  The  clas¬ 
sification  is  made  by  comparing  the  values  of  the  objective  function  Q  of  each  DoLP 
estimate  as  described  in  Section  2.6.1. 

Given  the  metal  and  dielectric  DoLP  priors,  the  first  update  equation  in  Eq.  (2.57) 
produces  two  DoLP  estimates  P.  To  compute  the  value  of  Q  for  each  material  class, 
these  DoLP  estimates  are  substituted  into  Eq.  (3.1),  which  is  obtained  from  Eqs. 
(2.56)  and  (2.58),  along  with  the  estimates  of  the  remaining  parameters  S0,  hc,  and 

ic: 


QH+1(P ,  S0,a,h)=  [Xit,u  In  fetor, u)  +  Xest.p  ln  (XptotJ. 


X 


EE 


dc(  y) 

ltc(y) 


K(y  ~xo) 


+  (Xest,u  +  Xest.p)  EE{; 

n  v  ' 


X  In  [K+\y  -  x0)]}  -  C  (Xest,u  +  Xest,p) 


dc(  y) 


l(y) 


^(y-xo)  (3.i) 


where 


=  \  [1  -  P"(xo)]  SJ(X») 

Xest,p  =  Pn(*o)SZ (xo)cos2  K(x0)  -  9C] 
x”;.1,,.  =  t1  -  S„”+1(x0) 

Xptotp  =  pptt,r(xo)^+1(xo)cos2  [an+1(x0)  -  ec\ 
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In  Eq.  (3.2),  the  subscripts  n  and  n  +  1  indicate  estimates  at  the  nth  (previous) 
and  (n  +  l)th  (current)  iterations,  P"(x0)  is  the  previous  DoLP  estimate  at  pixel  x0, 
■^priori X-o)  is  the  current  DoLP  estimate  at  x0  computed  using  the  metal  or  dielectric 
DoLP  prior,  and  C  is  the  total  number  of  polarization  channels.  The  DoLP  estimate 
(either  the  estimate  using  the  metal  DoLP  prior  or  the  estimate  using  the  dielectric 
DoLP  prior)  which  maximizes  Q  determines  whether  the  material  is  classified  as  a 
metal  or  a  dielectric  at  pixel  x0: 


material  at  xq= 


metal  Qd(x0)  ~  Qm(x0)  <  0 
dielectric  else 


(3.3) 


where  Qd(x0)  and  Qm.(x0 )  denote  the  Q  values  using  the  dielectric  and  the  metal 
DoLP  priors  at  pixel  x0,  respectively  [11], 


3.2  Adaptive  DoLP  Priors 

Section  2.6.2  mentions  that  natural  materials  weakly  polarize  scattered  light  at 
near-normal  incident  angles.  This  causes  Hyde’s  material-classification  algorithm  to 
perform  poorly  at  near-normal  collection  geometries  [11].  Figure  3.1  depicts  the  prob¬ 
ability  distributions  of  the  measured  DoLP  at  incident  angles  9i  =  50°  and  9i  =  24°. 
The  target  scene  consists  of  metal  and  dielectric  parts.  The  probability  distribution 
at  9i  =  50°  is  consistent  with  Hyde’s  DoLP  priors  shown  in  Fig.  2.5.  Contrast  this 
with  the  distribution  shown  in  Fig.  3.1(b).  The  range  of  DoLP  values  observed  in  the 
9i  =  24°  case  is  approximately  1/4  of  the  9i  =  50°  case,  which  demostrates  the  main 
limitation  of  Hyde’s  implementation.  The  threshold  between  the  metal  and  dielectric 
DoLP  values  at  the  9i  =  24°  geometry  becomes  so  low  that  Hyde’s  fixed  DoLP  prior 
models  are  not  able  to  accurately  classify  materials  in  this  scenario.  Hyde  suggests 
that  this  limitation  can  be  alleviated  by  adaptively  updating  the  DoLP  priors  as  more 
information  about  the  scene  becomes  available  [11],  The  following  sections  introduce 
three  novel  approaches  to  implement  adaptive  DoLP  priors,  termed  hereafter  higher- 
order  super-Gaussian  method,  Gaussian  method,  and  distribution-averaging  method. 
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Figure  3.1:  (a)  Distribution  of  measured  DoLP  at  6j  =  50°.  (b)  Distribution  of 

measured  DoLP  at  0*  =  24°. 


3.2.1  Higher- Order  Super- Gaussian  Method.  Hyde’s  pBRDF  study  finds 
that  the  DoLP  is  best  modeled  as  a  uniform  probability  density  between  0  and  1 
because  of  numerous  uncontrollable  factors  such  as  observation  geometry  and  surface 
state  [11].  In  order  to  derive  the  DoLP  update  equations  in  Eq.  (2.59),  the  uniform 
distribution  of  the  DoLP  prior  needs  to  be  modeled  with  a  differentiable  distribution 
for  which  a  continuous  super-Gaussian  distribution  in  Eq.  (2.54)  is  chosen  [11], 

In  the  higher-order  super-Gaussian  method,  the  sum  of  two  super-Gaussian 
distributions  is  fit  to  the  estimated  DoLP  distribution.  The  fit  function  is  of  the  form 


n(P)  =  cme"[a"l(p"bm)r  +  cde-[adiP~bd)r ,  (3.4) 

where  the  coefficients  cm  and  q  are  normalizing  factors,  am  and  determine  the 
widths  of  the  super-Gaussians,  bm  and  ba  determine  the  centers  of  the  super-Gaussians, 
and  m  is  the  order  [11,13].  The  subscripts  m  and  d  represent  a  metal  and  a  dielectric, 
respectively.  Figure  3.2(a)  shows  the  resulting  fit  of  Eq.  (3.4)  with  the  estimated 
DoLP  distribution  at  6i  =  24°.  Figure  3.2(b)  shows  the  resulting  adaptive  DoLP 
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Figure  3.2:  Higher-order  super-Gaussian  method  for  adaptive  DoLP  priors,  (a) 

Distribution  of  DoLP  measurements  and  the  fit  of  Eq.  (3.4)  when  m  =  10.  (b)  Metal 
and  dielectric  priors  formed  from  fit.  The  constants  in  Eq.  (3.4)  determined  by  the 
fit  are  cm  —  1,  am  —  64.3,  bm  =  0.0195  and  cd  —  1,  a^—  13.7,  bd  =  0.104. 

priors  obtained  from  the  fit  shown  in  Fig.  3.2(a).  The  obtained  coefficients  are  then 
substituted  into  the  update  equation  in  Eq.  (2.59)  to  compute  the  metal  and  dielectric 
DoLP  estimates.  The  pixel  is  then  classified  using  Eq.  (3.3). 

This  approach  is  analytically  sound,  but  has  a  major  drawback.  Figure  3.2  re¬ 
veals  that  the  higher-order  super-Gaussian  fitting  and  the  prior-formation  processes 
generate  a  threshold  which  is  significantly  lower  than  the  one  empirically  observed  in 
the  DoLP  distribution.  As  will  be  shown,  this  drawback  leads  to  significant  classifi¬ 
cation  errors. 


3.2.2  Gaussian  method.  Equation  (2.15)  states  that  the  DoLP  is  the  ratio 
of  the  linearly  polarized  irradiance  to  the  total  irradiance.  Thus,  assuming  that  shot 
noise  is  the  dominant  noise  source,  DoLP  measurements  observed  at  a  given  collection 
geometry  and  surface  state  should  follow  the  random  variable  distribution  formed  by 
dividing  two  Poisson  random  variables.  Unfortunately,  no  analytical  form  exists  for 
this  PDF.  Figure  3.3  shows  a  PDF  of  this  distribution  obtained  via  Monte  Carlo 
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Figure  3.3:  PDF  of  the  division  of  two  Poisson  random  variables  using  Monte 

Carlo  analysis.  The  plot  shows  the  PDF  of  DoLP  measurements  of  a  target  whose 
true  DoLP  is  0.5  (SNR  =  \/l000). 


analysis.  The  figure  shows  the  distribution  of  DoLP  measurements  of  a  target  whose 
true  DoLP  is  0.5  (SNR  =  \/l000).  Note  that  the  DoLP  PDF  is  very  Gaussian-like. 
Based  on  this  finding,  it  seems  much  more  reasonable  to  fit  the  sum  of  two  Gaussian 
functions  to  the  estimated  DoLP  distribution  than  it  does  super-Gaussians.  Thus,  in 
the  Gaussian  method, 


n  (P)  =  cme-[am(p-fem)]2  +  cde-[adiP~bd)]2  (3.5) 

is  fit  to  the  measured  DoLP  PDF. 

Fig.  3.4(a)  shows  the  result  of  this  fit  for  DoLP  measurements  at  9j  =  24°.  The 
figure  shows  that  the  Gaussian  fit  models  DoLP  measurements  more  accurately  than 
super-Gaussian  fit  depicted  in  Fig.  3.2(a).  This  Gaussian  fit  provides  the  threshold 
between  a  metal  and  a  dielectric  from  which  the  DoLP  priors  can  be  built.  As 
mentioned  in  the  previous  section,  the  super-Gaussian  distribution  is  chosen  to  model 
the  uniform  distribution  of  the  overall  DoLP.  In  order  to  compute  the  coefficients 
in  Eq.  (3.4)  necessary  to  model  the  updated  DoLP  priors  using  the  super-Gaussian 
distribution,  the  threshold  between  the  two  materials  are  determined  by  finding  the 
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Figure  3.4:  Gaussian  method  for  adaptive  DoLP  priors,  (a)  Distribution  of  DoLP 
measurements  and  the  fit  of  Eq.  (3.5).  (b)  Metal  and  dielectric  priors  formed  from  fit. 
The  constants  in  Eq.  (3.5)  determined  by  the  fit  are  cm  =  1 ,  am  =  34.5,  bm  =  0.0315 
and  q  =  1 ,  a,d  =  2.32,  bd  =  0.532. 

minimum  between  the  two  maxima  of  the  Gaussian  fit  (marked  in  Fig.  3.4(a)).  The 
center  of  each  prior  is  then  obtained  by  computing  the  midpoint  between  the  threshold 
and  the  DoLP  values  of  0  and  1.  Given  the  threshold  and  the  center  of  a  prior,  the 
coefficient  a  is  computed  using  Eq.  (2.54),  i.e. , 


[-  lnn(Ptfc)]1/i 

(Pth  -  b ) 


(3.6) 


where  Pth  is  the  DoLP  value  at  the  threshold  and  fl(P^)  is  a  user-chosen  quantity.  In 
this  analysis,  its  value  is  chosen  to  be  0.01.  Figure  3.4(b)  shows  the  updated  DoLP 
priors  with  the  coefficients  computed  as  described. 


3.2.3  Distribution- Averaging  Method.  The  final  method  investigated  for 
adaptively  updating  the  DoLP  priors  is  to  approximate  the  threshold  by  Ending  the 
mean  of  the  measured  DoLP  distribution.  It  has  been  found  empirically  that  the 
mean  of  the  distribution  lies  approximately  at  the  intersection  (or  threshold)  of  the 
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Figure  3.5:  Distribution-averaging  method  for  adaptive  DoLP  priors,  (a)  Distribu¬ 
tion  of  DoLP  measurements  and  the  calculated  mean,  (b)  Updated  DoLP  priors.  The 
constants  in  Eq.  (3.4)  determined  by  the  mean  are  cm  =  1,  am  =  45.4,  bm  =  0.0257 
and  Cd  =  1,  ad  =  2.46,  bd  =  0.109. 

two  “humps”  of  the  DoLP  distribution.  In  order  to  update  the  DoLP  priors,  the 
coefficients  of  the  super-Gaussian  distribution  in  Eq.  (3.4)  are  computed  from  the 
threshold  and  the  DoLP  values  of  0  and  1,  using  the  same  process  as  in  the  Gaussian 
method.  Figure  3.5(b)  shows  the  resulting  metal  and  dielectric  DoLP  priors. 


3.3  Algorithm  Execution 

Figure  3.6  outlines  a  flow  of  the  proposed  material-classification  algorithm  at 
the  ith  iteration.  The  number  of  iterations  before  forming  the  priors  and  classifying 
the  object  is  determined  by  comparing  the  DoLP  estimates  (obtained  using  standard 
LeMaster  and  Cain  deconvolution)  to  the  DoLP  obtained  from  turbulence- free  polari- 
metric  images.  The  root  mean  squared  error  (RMSE)  of  DoLP  estimates  in  Fig.  3.7(a) 
indicates  that  the  DoLP  estimate  is  optimized  at  about  50  iterations.  Thus,  on  the 
51st  iteration,  the  DoLP  priors  are  updated  using  one  of  three  techniques  developed 
in  this  chapter  and  the  binary  classification  of  the  objects  in  the  scene  commences. 
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Figure  3.6:  Flowchart  of  the  proposed  material-classification  algorithm  augmented 
with  adaptive  DoLP  priors  at  the  ith  iteration. 

The  stopping  criterion  for  the  algorithm  is  decided  based  on  the  Sq  estimation 
results.  Shown  in  Fig.  3.7(b)  is  the  RMSE  of  the  So  estimation  where  the  irradiance 
of  the  object  measured  without  atmospheric  turbulence  is  considered  the  true  object 
irradiance.  Note  that  the  turbulence-free  images  do  contain  static  aberrations  due 
to  the  imperfect  imaging  system.  This  explains  why  the  Sq  RMSE  increases  with 
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Figure  3.7:  (a)  RMSE  of  DoLP  to  determine  the  iteration  number  to  update 

DoLP  priors,  (b)  RMSE  of  So  t°  determine  the  stopping  criterion  for  the  material- 
classification  algorithm. 


the  number  of  iterations  before  converging.  Based  on  the  observation  that  the  So 
estimation  stagnates  after  300  iterations,  the  entire  algorithm  is  stopped  after  300 
iterations. 
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IV.  Results  and  Analisys 


This  chapter  presents  results  of  the  material-classification  algorithm  introduced 
in  Chapter  3.  Section  4.1  describes  the  experimental  design  and  procedures. 
Sections  4.2  and  4.3  discuss  classification  results  of  the  algorithm  using  fixed  and 
adaptive  DoLP  priors.  A  spatial  error  analysis  of  the  results  is  also  presented  in 
Section  4.4. 

4-1  Experimental  Setup 

The  images  analyzed  in  this  research  were  collected  by  Hyde  using  a  Stokes 
polarimeter  in  the  Optical  Turbulence  Estimation,  Compensation,  and  Simulation 
laboratory  at  the  Air  Force  Institute  of  Technology.  Figure  4.1  shows  the  instrument, 
which  consists  of  three  main  sections.  The  source/sample  section  consists  of  two 
optical  rotation  stages.  The  lower  rotation  stage  holds  the  sample  while  the  upper 
rotation  stage,  which  is  placed  on  the  top  of  the  lower  rotation  stage,  holds  the  1550 
nm  light  emitting  diode  source.  These  two  stages  are  rotated  in  combination  to  yield 
any  9i  +  6r  >  48°  where  6i  and  9r  denote  the  incident  and  reflected  angles,  respec¬ 
tively.  Aluminum  and  steel  samples  are  partially  painted  to  simulate  a  target  scene 


Polarizer/Camera 


Turbulence  Simulator 
Phase  Wheel 


1550  nm  LED  Jkl. 
fiaSy  Sample  Holder 


Samples 


Figure  4.1:  Photograph  of  the  Stokes  polarimeter  used  to  collect  the  polarimetric 
imagery  in  this  experiment  [11]. 
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consisting  of  metal  and  dielectric  parts.  The  atmospheric  turbulence  simulator  section 
of  the  instrument  models  various  atmospheric  turbulence  conditions  by  changing  the 
location  of  a  phase  wheel.  Atmospheric  turbulence  conditions  used  in  this  research 
are  D /r0  ~  7.9,  10.5,  and  12.9.  The  last  section  of  the  instrument  is  the  polarimetric 
imager  which  includes  an  imaging  lens  (pupil  size  of  25  mm),  a  rotating  linear  polar¬ 
izer,  and  a  digital  camera  whose  array  size  is  250  x  316  with  a  pixel  pitch  of  30  fim. 
For  more  details  about  the  experimental  setup,  the  reader  is  referred  to  [11]. 

In  order  to  implement  the  proposed  material-classification  algorithm,  the  Stokes 
parameters  S0,  Si,  and  S2  are  extracted  from  four  polarimetric  images  ( 9C  =  0°,  45°, 
90°,  135°).  Dark- frame  images  were  collected  and  subtracted  from  the  polarimetric 
images  before  computing  the  Stokes  parameters.  The  original  images  of  size  250  x 
316  are  padded  to  size  512  x  512  and  then  spatially  windowed  using  a  Tukey  window 
to  prevent  “ringing”  caused  by  the  fast  Fourier  transforms  used  in  the  algorithm. 

The  algorithm  starts  with  choosing  initial  estimates  of  the  sought  parameters — 
P,  Sq,  and  the  PSF  hc.  The  initial  estimates  of  P  and  So  are  512  x  512  arrays  of 
1/2  and  1,  respectively.  The  initial  estimate  of  hc  is  formed  from  a  randomly  drawn 
atmospheric  phase  screen  ( D/ro  =  10)  with  piston  and  tilt  removed. 

4-2  Classification  Results — Fixed  DoLP  Priors 

Figure  4.2  shows  results  of  the  material-classification  algorithm  with  Hyde’s 
DoLP  priors  for  the  painted  aluminum  target  after  300  iterations.  The  atmospheric 
turbulence  condition  and  the  collection  geometry  are  D/ro  ~  7.9  and  =  Qr  =  50°, 
respectively.  Figures  4.2(a),  (d),  and  (g)  show  S0  images  of  no-turbulence,  turbulence- 
degraded,  and  estimated  data,  respectively.  The  experimental  P  results  are  shown  in 
Figs.  4.2(b),  (e),  and  (h).  The  threshold  of  P  =  0.3  is  chosen  to  discriminate  between 
a  metal  (P  <  0.3)  and  a  dielectric  (P  >  0.3)  for  the  no-turbulence  and  turbulence- 
degraded  results.  This  threshold  value  is  the  intersection  of  the  two  DoLP  priors  as 
shown  in  Fig.  2.5.  The  corresponding  material-classification  results  are  shown  in  Figs. 
4.2(c),  (f),  and  (i).  As  indicated  by  Hyde’s  previous  research,  the  estimated  classifi- 
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Figure  4.2:  Results  of  Hyde’s  material-classification  algorithm  for  the  painted  alu¬ 
minum  target  {D/tq  m  7.9  and  6t  —  9r  —  50°).  For  the  no-turbulence  and  the 
turbulence-degraded  measurements,  P  >  0.3  is  classified  as  a  dielectric  and  P  <  0.3 
as  a  metal.  The  estimated  classification  result  is  based  on  the  metal  and  dielectric 
DoLP  priors  developed  by  Hyde  and  reviewed  in  Section  2.6.1.  In  the  classification 
results,  the  blue  color  implies  dielectric  and  the  red  color  implies  metal.  The  final 
estimates  are  obtained  after  300  iterations. 


cation  result  in  Fig.  4.2 (i)  reveals  that  the  algorithm  not  only  accurately  classifies  the 
object  but  also  recovers  details  of  the  target  scene  lost  by  atmospheric  turbulence, 
evident  in  the  estimated  P  result. 

Figure  4.3,  on  the  other  hand,  illustrates  the  limitation  of  the  algorithm  using 
fixed  DoLP  priors  for  near- normal  collection  geometries  (9t  =  9r  =  24°  in  this  case). 
The  threshold  of  P  =  0.08  is  chosen  from  the  histogram  of  the  DoLP  measurements 
in  Fig.  3.1(b)  to  discriminate  between  a  metal  (P  <  0.08)  and  a  dielectric  (P  >  0.08) 
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Figure  4.3:  Results  of  Hyde’s  material-classification  algorithm  for  the  painted  alu¬ 
minum  target  ( D/r0  m  7.9  and  9t  =  9r  =  24°).  For  the  no-turbulence  and  the 
turbulence-degraded  measurements,  P  >  0.08  is  classified  as  a  dielectric  and  P  <  0.08 
as  a  metal.  The  estimated  classification  result  is  based  on  the  metal  and  dielectric 
DoLP  priors  developed  by  Hyde  and  reviewed  in  Section  2.6.1.  In  the  classification 
results,  the  blue  color  implies  dielectric  and  the  red  color  implies  metal.  The  final 
estimates  are  obtained  after  300  iterations. 


for  the  no-turbulence  and  turbulence-degraded  data.  To  demonstrate  the  limitation 
of  the  algorithm  in  this  near-normal  geometry,  the  DoLP  priors  developed  by  Hyde 
and  reviewed  in  Section  2.6.1  are  incorporated  into  the  algorithm  for  material  classi¬ 
fication.  The  estimated  classification  result  in  Fig.  4.3(i)  confirms  the  predicted  poor 
performance  of  the  algorithm  for  near-normal  collection  geometries. 
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4-3  Classification  Results — Adaptive  DoLP  Priors 

In  order  to  alleviate  the  near-normal  collection  geometry  limitation  of  the  fixed- 
prior  material-classification  algorithm,  the  DoLP  priors  should  be  updated  as  more 
information  about  the  target  scene  becomes  available.  The  following  sections  present 
classification  results  of  the  three  DoLP  prior  update  methods  developed  in  Chapter  III. 

4-3.1  Higher-Order  Super- Gaussian  Method.  Figure  4.4  shows  results  of 
the  material-classification  using  adaptive  DoLP  priors  updated  by  the  higher-order 
super-Gaussian  method  discussed  in  Section  3.2.1.  The  same  images  of  the  painted 
aluminum  target  (D/r o  ~  7.9,  dj  =  9r  =  24°)  in  Fig.  4.3  are  used  for  the  proposed 
method.  In  order  to  classify  the  no-turbulence  and  the  turbulence-degraded  measure¬ 
ments,  the  threshold  value  is  computed  by  fitting  the  distribution  of  each  DoLP  mea¬ 
surement  with  the  higher-order  super-Gaussian  function  as  described  in  Section  3.2.1. 
The  resulting  thresholds  of  P  =  0.0618  and  P  =  0.0464  are  then  applied  to  classify 
the  no-turbulence  and  the  turbulence-degraded  data,  respectively. 

The  estimated  classification  result  using  the  updated  DoLP  priors  is  shown 
in  Fig.  4.4(i).  Note  that  the  figure  depicts  only  a  few  pixels  correctly  classified  as 
metal.  Figure  4.4(i)  indicates  that  the  algorithm  with  DoLP  priors  updated  by  the 
higher-order  super-Gaussian  method  performs  as  poorly  as  the  algorithm  with  fixed 
DoLP  priors  shown  in  Fig.  4.3(i).  As  mentioned  in  Section  3.2.1,  the  higher-order 
super-Gaussian  method  yields  a  lower  threshold  than  the  one  empirically  observed  in 
the  measured  DoLP  distribution.  This  causes  the  incorrect  dielectric  over  classifica¬ 
tion.  Although  the  estimated  classification  result  of  the  higher-order  super-Gaussian 
method  is  poor,  the  blind-deconvolution  algorithm  still  produces  Sq  and  P  estimates 
with  recovered  details  lost  in  the  turbulence-degraded  images. 

4-3.2  Gaussian  Method.  Figure  4.5  shows  results  of  the  material-classification 
algorithm  using  DoLP  priors  adaptively  updated  by  the  Gaussian  method  discussed 
in  Section  3.2.2.  The  images  of  the  target  and  the  layout  of  the  figure  are  identical  to 
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Figure  4.4:  Results  of  the  material-classification  algorithm  using  adaptive  DoLP 

priors  derived  from  the  higher-order  super- Gaussian  method  for  the  painted  aluminum 
target  ( 9j  =  9r  =  24°  and  D/r0  ss  7.9).  The  thresholds  for  classifying  the  no- 
turbulcnce  and  the  turbulence-degraded  measurements  are  P  =  0.0618  and  P  = 
0.0464,  respectively.  The  estimated  classification  result  is  based  on  the  metal  and 
dielectric  DoLP  priors  adaptively  updated  by  the  higher-order  super-Gaussian  method 
discussed  in  Section  3.2.1.  In  the  classification  results,  the  blue  color  implies  dielectric 
and  the  red  color  implies  metal.  The  final  estimates  are  obtained  after  300  iterations. 


the  previous  results.  Similar  to  the  higher-order  super-Gaussian  method,  to  classify 
the  no-turbulence  and  the  turbulence-degraded  images,  the  threshold  is  computed 
by  fitting  the  distribution  of  the  DoLP  measurements  with  the  Gaussian  function  as 
described  in  Section  3.2.2.  The  resulting  thresholds  are  P  =  0.0932  and  P  =  0.0747 
for  the  no-turbulence  and  turbulence-degraded  data,  respectively. 
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Figure  4.5:  Results  of  the  material-classification  algorithm  using  adaptive  DoLP  pri¬ 
ors  derived  from  the  Gaussian  method  for  the  painted  aluminum  target  (#*  =  6r  =  24° 
and  D/r0  m  7.9).  The  thresholds  for  classifying  the  no-turbulence  and  the  turbulence- 
degraded  measurements  are  P  =  0.0932  and  P  =  0.0747,  respectively.  The  estimated 
classification  result  is  based  on  the  metal  and  dielectric  DoLP  priors  adaptively  up¬ 
dated  by  the  Gaussian  method  discussed  in  Section  3.2.2.  In  the  classification  results, 
the  blue  color  implies  dielectric  and  the  red  color  implies  metal.  The  final  estimates 
are  obtained  after  300  iterations. 


Figure  4.5(i)  shows  that  the  algorithm  using  adaptive  priors  updated  by  the 
Gaussian  method  accurately  classifies  the  target  scene.  Shown  in  Fig.  4.5(h)  is  the 
estimated  P  result.  Note  the  clear  contrast  between  the  two  material  types  which 
comprise  the  target  scene.  This  result  is  a  clear  validation  of  this  approach. 

The  estimated  So  result  is  similar  to  the  result  of  the  previous  method.  Although 
the  S'o  result  seems  to  recover  features  lost  in  the  turbulence-degraded  images,  it  is  dif¬ 
ficult  to  objectively  quantify  this  observation  because  of  the  aberrated  no-turbulence 
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Figure  4.6:  Results  of  the  material-classification  algorithm  using  adaptive  DoLP 

priors  derived  from  the  distribution-averaging  method  for  the  painted  aluminum  tar¬ 
get  [6i  =  9r  =  24°  and  D/r0  ~  7.9).  The  thresholds  for  classifying  the  no-turbulence 
and  the  turbulence-degraded  measurements  are  P  =  0.1223  and  P  =  0.0651,  respec¬ 
tively.  The  estimated  classification  result  is  based  on  the  metal  and  dielectric  DoLP 
priors  adaptively  updated  by  the  distribution-averaging  method  discussed  in  Section 
3.2.3.  In  the  classification  results,  the  blue  color  implies  dielectric  and  the  red  color 
implies  metal.  The  final  estimates  are  obtained  after  300  iterations. 


images.  As  Section  3.3  briefly  mentions,  the  no-turbulence  images  which  are  consid¬ 
ered  “truth”  in  this  research  for  the  basis  of  comparison  do  contain  static  imaging 
aberrations. 


4-3.3  Distribution- Averaging  Method.  Figure  4.6  shows  results  of  the  material- 

classification  algorithm  using  the  distribution-averaging  method  developed  in  Sec¬ 
tion  3.2.3.  The  images  of  the  target  and  the  layout  of  the  figure  are  identical  to  those 
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of  the  previous  results.  To  classify  the  no-turbulence  and  the  turbulence-degraded 
data,  the  threshold  is  computed  by  finding  the  mean  of  the  distribution  of  the  DoLP 
measurements.  The  corresponding  thresholds  are  P  =  0.1223  and  P  =  0.0651  for  the 
no-turbulence  and  the  turbulence-degraded  data,  respectively. 

The  estimated  classification  result  in  Fig.  4.6(i)  shows  that  DoLP  priors  updated 
by  the  distribution-averaging  method  significantly  improve  the  performance  of  the 
fixed-threshold  based  algorithm  at  the  tested  near-normal  collection  geometry.  Per¬ 
formance  is  similar  to  the  previously  analyzed  Gaussian  method  with  some  instances 
of  misclassihcation  occurring  along  the  metal/dielectric  boundaries.  This  may  be  an 
acceptable  trade-off  considering  the  computational  simplicity  of  this  technique. 

Like  the  Gaussian  method  result,  the  estimated  P  result  in  Fig.  4.6(h)  displays 
a  clear  contrast  between  the  two  material  classes  with  features  recovered  after  being 
lost  in  the  turbulence-degraded  images.  The  estimated  S0  result  does  not  display  any 
noticeable  differences  from  the  S'o  estimation  results  of  the  previous  methods.  Over¬ 
all,  the  distribution-averaging  method  yields  good  results  for  such  a  computationally 
simplistic  technique. 

4-4  Spatial  Analysis  of  the  Classification  Results 

Table  4.1  summarizes  the  material-classification  results  of  the  three  proposed 
methods  for  the  painted  aluminum  target  at  the  near-normal  collection  geometry 
tested  in  this  experiment.  The  table  reports  the  results  for  three  atmospheric  turbu¬ 
lence  conditions  ( D/r0  ps  7.9, 10.5,  and  12.9).  The  first  column  indicates  the  method 
applied  to  update  the  DoLP  priors.  The  second  column  indicates  the  image  being 
analyzed,  either  turbulence-degraded  (Turb.)  or  estimated  (Est.)  images.  The  last 
two  columns  contain  the  RMSE  and  the  correlation  coefficient  of  the  classification 
results. 
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Table  4.1: 
9r  =  24° 


Material-classification  results  for  the  painted  aluminum  target  at 


Adaptive  Method 

Image 

D/r  o 

RMSE 

Corr.  Coeff. 

Turb. 

7.9 

0.575 

0.336 

Higher- Order 
Super-Gaussian 

Est. 

Turb. 

7.9 

10.5 

0.730 

0.581 

0.0337 

0.306 

Method 

Est. 

10.5 

0.733 

0.0080 

Turb. 

12.9 

0.595 

0.286 

Est. 

12.9 

0.734 

0.0095 

Turb. 

7.9 

0.396 

0.562 

Est. 

7.9 

0.391 

0.573 

Gaussian 

Turb. 

10.5 

0.399 

0.555 

Method 

Est. 

10.5 

0.396 

0.563 

Turb. 

12.9 

0.405 

0.542 

Est. 

12.9 

0.409 

0.534 

Turb. 

7.9 

0.369 

0.536 

Distribution- 

Est. 

7.9 

0.376 

0.529 

Turb. 

10.5 

0.383 

0.514 

Averaging 

Method 

Est. 

10.5 

0.388 

0.509 

Turb. 

12.9 

0.398 

0.491 

Est. 

12.9 

0.393 

0.498 

The  nearly-zero  correlation  coefficients  of  the  higher-order  super-Gaussian  method 
in  Table  4.1  conclusively  show  that  this  technique  fails  to  accurately  classify  the  target 
for  the  reason  discussed  previously. 

The  Gaussian  method,  whose  computed  threshold  is  the  best  approximation  to 
empirical  data  among  the  three  methods,  shows  slightly  improved  performance  over 
the  classification  using  turbulence-degraded  images.  These  improvements,  however, 
are  quantitatively  too  minimal  to  validate  the  proposed  algorithm.  Thus,  in  order  to 
evaluate  the  validity  of  the  algorithm  properly,  it  needs  to  be  understood  that  the 
no-turbulence  images  (“truth”)  contain  static  aberrations  as  previously  mentioned. 
These  aberrations  penalize  the  estimated  classification  result  when  computing  the 
RMSE  and  the  correlation  coefficient  because  the  algorithm  corrects  these  aberra- 
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tions  as  qualitatively  shown  in  Figs.  4.5(h)  and  (i).  Despite  marginal  quantitative 
improvements,  the  estimated  P  and  classification  results  shown  in  Figs.  4.5(h)  and  (i) 
clearly  depict  that  the  Gaussian  method  classifies  the  target  with  greater  accuracy 
than  the  turbulence-degraded  images.  Further  analysis  using  no-turbulence  images 
whose  static  aberrations  are  removed  may  provide  more  relevant  image  metrics  for 
quantitative  study  of  the  classification  results. 

The  metric  results  of  the  distribution-averaging  method  are  similar  to  those  of 
the  Gaussian  method  in  that  improvements  in  performance  increase  as  atmospheric 
turbulence  becomes  weaker  ( D/r0  decreases).  Unlike  the  Gaussian  update  tech¬ 
nique,  however,  the  distribution-averaging  method  yields  classification  results  which 
are  quantitatively  worse  than  the  turbulence-degraded  results.  For  the  same  rea¬ 
son  discussed  above,  it  can  still  be  concluded  that  the  distribution-averaging  method 
improves  the  classification  over  the  turbulence-degraded  result  based  on  Figs.  4.6(h) 
and  (i).  Note  that  Table  4.1  indicates  that  the  Gaussian  method  provides  results  some¬ 
what  better  than  the  distribution-averaging  method.  Considering  the  computational 
simplicity  of  the  latter  update  technique,  relatively  small  differences  in  performance 
may  be  an  acceptable  trade-off  as  suggested  in  Section  4.3.3. 
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V.  Conclusions  and  Future  Work 


5. 1  Conclusions 

In  this  thesis,  a  binary  material-classification  algorithm  augmented  with  adap¬ 
tive  DoLP  priors  is  developed  for  improved  performance  at  near-normal  collection 
geometries.  The  algorithm  is  based  on  Hyde’s  material-classification  technique  which 
utilized  fixed  DoLP  priors.  Hyde’s  algorithm  used  a  variant  of  the  multichannel  blind- 
deconvolution  technique  developed  by  LeMaster  and  Cain  to  estimate  So,  P ,  and  hc 
from  polarimetric  imagery  degraded  by  atmospheric  turbulence.  To  discriminate  be¬ 
tween  metals  and  dielectrics,  DoLP  estimates  from  the  blind-deconvolution  results 
of  each  DoLP  prior  are  used  to  compute  the  corresponding  values  of  the  objective 
function  Q.  Whichever  DoLP  prior  maximizes  0  determines  the  classification  of  the 
target  at  a  given  pixel  location.  A  key  limitation  of  this  algorithm  arises  from  the 
polarization  behavior  of  light  scattered  from  natural  materials  at  near-normal  collec¬ 
tion  geometries.  The  previously  developed  material-classification  technique  relied  on 
fixed-threshold  DoLP  priors  which  performed  poorly  for  near-normal  collection  ge¬ 
ometries.  To  alleviate  this  limitation  in  collection  geometries,  this  research  developed 
methods  to  adaptively  update  the  DoLP  priors  as  more  information  about  the  scene 
became  available. 

Three  methods  are  investigated  for  performing  this  update.  The  first  method, 
called  the  higher-order  super- Gaussian  method,  involves  fitting  the  sum  of  two  super- 
Gaussian  functions  to  the  measured  distribution  of  DoLP  estimates.  This  method  is 
analytically  sound,  but  suffers  because  the  resulting  fit  poorly  predicts  the  location  of 
the  threshold  empirically  evident  in  the  estimated  DoLP  distributions.  This  results  in 
misclassification  of  metal  as  dielectric  at  numerous  pixels.  The  second  method,  called 
the  Gaussian  method,  fits  the  distribution  of  DoLP  estimates  with  the  sum  of  two 
Gaussian  functions.  This  approach  is  based  on  the  observation  that  DoLP  measure¬ 
ments  (at  a  given  surface  state  and  collection  geometry)  are  equivalent  to  a  Poisson 
random  variable  divided  by  another  Poisson  random  variable.  This  distribution  has 
no  analytical  form  but  is  Gaussian-like.  It  is  clearly  shown  that  the  Gaussian  method 
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computes  the  most  accurate  threshold.  The  last  technique  developed  and  analyzed  in 
this  thesis  is  the  distribution-averaging  method.  This  method  is  based  on  an  empirical 
observation  that  the  mean  of  the  DoLP  distribution  approximately  lies  at  the  true 
threshold  value.  This  simple  approach  produces  classification  results  which  are  signif¬ 
icantly  improved  over  the  existing  fixed-prior  method.  Both  RMSE  and  correlation 
coefficient  calculations  of  experimental  data  confirm  that  the  DoLP  priors  updated 
by  the  Gaussian  method  and  the  distribution-averaging  method  produce  accurate 
classification  results  in  the  tested  atmospheric  turbulence  conditions. 

5. 2  Future  Work 

5.2.1  Testing  Collection  Geometry  Region  of  Validity.  The  first  area  of  pos¬ 
sible  future  work  is  to  perform  experiments  on  the  proposed  methods  with  collection 
geometries  9%,  6r  <  24°.  Even  though  the  adaptive  approach  developed  in  this  research 
significantly  extends  the  region  of  applicability  of  the  algorithm  beyond  the  existing 
method,  the  algorithm  will  eventually  fail  because  the  polarization  difference  between 
dielectrics  and  metals  will  become  too  small  to  be  discerned.  The  collection  geometry 
range  of  validity  for  the  proposed  methods  has  not  been  confirmed  empirically.  Ex¬ 
periments  with  collection  geometries  9^,  9r  <  24°  will  allow  one  to  fully  understand 
the  capabilities  of  the  proposed  techniques. 

5.2.2  Field  Testing.  Next  the  proposed  algorithm  should  be  applied  to  a 
real-world  scenario  by  conducting  field  experiments.  A  compact  polarimeter  would 
be  needed  for  this  purpose.  Figure  5.1  depicts  a  polarization  camera  suitable  for  such 
field  testing,  which  is  also  referred  to  in  Hyde’s  research.  This  camera  is  a  commercial 
division  of  amplitude  polarimeter  with  three  polarization  channels  configured  with  0°, 
45°,  and  90°  linear  polarizers  to  measure  the  Stokes  parameters  S0,  Si,  and  S2. 

5.2.3  Comparison  of  Material-Classification  Results.  There  exist  other 
polarization-based  material-classification  techniques,  such  as  Wolff’s  polarization  Fres¬ 
nel  ratio  and  Tominaga  and  Kimachi’s  degree  of  polarization  (DoP)  map  methods. 
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Figure  5.1:  A  three-channel  polarization  camera  developed  by  FluxData,  Inc.  [5]. 


Wolff  utilizes  the  Fresnel  reflection  coefficients  to  estimate  the  polarization  Fresnel 
ratio  whose  value  is  used  to  discriminate  between  metals  and  dielectrics  [27].  Torni- 
naga  and  Kimachi  estimate  the  DoP  from  the  irradiance  through  a  polarizer  and  use 
the  curvature  of  the  DoP  map  to  discriminate  between  the  two  material  classes  [24], 
By  comparing  the  results  from  these  techniques,  the  relative  accuracy /validity  of  the 
proposed  methods  can  be  assessed. 


5.2.4  Use  of  So  Estimate  to  Aid  Classification.  One  thing  not  used  in 
this  research  was  the  information  contained  in  the  Sq  estimate.  The  rather  obvious 
observation  that  metals  tend  to  scatter  light  more  strongly  than  dielectrics  could  be 
used  to  aid  the  polarization-based  classification  and  extend  the  collection  geometry 
range  of  validity  of  the  algorithm  even  further. 


5.2.5  DoLP  Derivative  Classification.  While  DoLP  varies  significantly  ver¬ 
sus  incidence/observation  angle  and  surface  state,  the  rate  of  change  of  the  DoLP 
remains  fairly  consistent  over  a  wide  range  of  conditions.  Furthermore,  the  deriva¬ 
tives  of  the  DoLPs  of  metals  and  dielectrics  are  very  distinct  with  metals  having  a 
DoLP  derivative  of  approximately  zero  and  dielectrics  having  a  large  DoLP  derivative 
(especially  versus  incidence/observation  angle).  This  trend  holds  even  for  very  near¬ 
normal  collection  geometries.  For  this  approach  to  be  implemented,  the  deconvolution 
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algorithm  must  be  reworked  to  estimate  the  derivative  of  the  DoLP  or  measurements 
at  two  slightly  different  observation  angles  need  to  be  collected. 

5.2.6  Diffraction- Limited  Truth  Images.  The  quantitative  analysis  of  classi¬ 
fication  results  shows  that  the  no-turbulence  ( “truth” )  images  are  of  poor  quality  such 
that  improvements  in  estimates  brought  about  by  the  proposed  algorithm  are  not  ev¬ 
ident  in  the  objective  measures.  Using  “diffraction-limited”  no-turbulence  images  as 
truths  may  provide  more  relevant  quantitative  measures.  This  can  be  accomplished 
in  two  ways.  The  first,  of  course,  is  to  collect  new  polarimetric  images.  This  is 
now  possible  with  the  new  polarimetric  imager  in  the  Optical  Turbulence  Estimation, 
Compensation,  and  Simulation  laboratory.  This  instrument  unfortunately  was  un¬ 
available  for  this  thesis.  The  second  is  to  deconvolve  the  no-turbulence  images.  This 
was  not  done  in  this  analysis  to  remain  consistent  with  the  algorithm  comparison, 
performed  by  Hyde  in  [11], 
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14.  ABSTRACT 

This  research  develops  an  enhanced  material-classification  algorithm  to  discriminate  between  metals  and  dielectrics  using 
passive  polarimetric  imagery  degraded  by  atmospheric  turbulence.  To  improve  the  performance  of  the  existing  technique 
for  near-normal  collection  geometries,  the  proposed  algorithm  adaptively  updates  the  degree  of  linear  polarization 
(DoLP)  priors  as  more  information  becomes  available  about  the  scene.  Three  adaptive  approaches  are  presented.  The 
higher-order  super-Gaussian  method  fits  the  distribution  of  DoLP  estimates  with  a  sum  of  two  super-Gaussian  functions 
to  update  the  priors.  The  Gaussian  method  computes  the  classification  threshold  value,  from  which  the  priors  are 
updated,  by  fitting  the  distribution  of  DoLP  estimates  with  a  sum  of  two  Gaussian  functions.  Lastly,  the 
distribution-averaging  method  approximates  the  threshold  value  by  finding  the  mean  of  the  DoLP  distribution.  The 
experimental  results  confirm  that  the  new  adaptive  method  significantly  extends  the  collection  geometry  range  of  validity 
for  the  existing  technique. _ 
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